home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 5 / MacMania 5.toast / / Tools&Utilities / Plotfoil 3.2 / naca-1.0 / naca6.f < prev    next >
Text File  |  1995-09-13  |  10KB  |  427 lines
  1. c
  2. c-----------------------------------------------------------------------------
  3. c
  4. c  Naca6.f -- contains routines to generate the NACA 1 and 6 series
  5. c             foils.
  6. c
  7. c  Written By: S.E.Norris
  8. c
  9. c  norris@cfd.mech.unsw.edu.au
  10. c
  11. c  $RCSfile: naca6.f,v $
  12. c  $Author: norris $
  13. c  $Revision: 1.5 $
  14. c  $Date: 1995/08/31 11:05:52 $
  15. c
  16. c  $Log: naca6.f,v $
  17. c  Revision 1.5  1995/08/31  11:05:52  norris
  18. c  *** empty log message ***
  19. c
  20. c  Revision 1.4  1995/08/31  07:30:19  norris
  21. c  Cleaned up code.
  22. c
  23. c  Revision 1.3  1995/08/31  03:35:01  norris
  24. c  Changed to use the directory.h header file that holds the directory name.
  25. c
  26. c  Revision 1.2  1995/08/31  02:03:05  norris
  27. c  Added the use of the lnblnk function to generate pathnames.
  28. c
  29. c  Revision 1.1  1995/08/30  11:07:20  norris
  30. c  Initial revision
  31. c
  32. c-----------------------------------------------------------------------------
  33. c
  34.       SUBROUTINE Naca1( x,npl,npmx,naca,inaca,scle,t_n )
  35. c
  36.       IMPLICIT none
  37.       INTEGER npl,npmx,inaca
  38.       REAL x(3,npmx),scle
  39.       CHARACTER naca*(*)
  40.       LOGICAL t_n
  41. c
  42. c  This routine calculates the ordinates of a NACA 1 series foil.
  43. c
  44.       INTEGER nbmx
  45.       PARAMETER (nbmx=30)
  46.       REAL pi
  47.       PARAMETER (pi=3.14159265358979)
  48. c
  49.       INTEGER npt,npf,nb,nbm
  50.       INTEGER i,j
  51.       REAL ndata(3)
  52.       REAL ab,bb,rb,xb(nbmx),yb(nbmx),y2b(nbmx)
  53.       REAL xi,xm
  54.       REAL dx,dy,dydx,t,y
  55.       CHARACTER filnme*6
  56.       REAL Nacat
  57.       EXTERNAL Nacat
  58. c
  59.       DATA  filnme / '16-015' /
  60. c
  61.       npt  = npl+1
  62.       npf  = npl/2
  63. c
  64. c     Calculate the parameters describing the foil, from the given NACA
  65. c     number. Then calculate the cubic spline for the camber.
  66. c
  67.       call Naca1b( naca,inaca,ndata )
  68.       call GetSpl( xb,yb,y2b,nb,nbm,nbmx,ab,bb,rb,filnme,t_n,ndata(3) )
  69.       if (.not. t_n) return
  70. c
  71. c     Set the x values according to a cosine distribution.
  72. c
  73.       call Cosdn( x,npl,npmx )
  74. c
  75. c     Calculate 'y' values for each value of 'x'.
  76. c
  77.       x(2,1)     = 0.0
  78.       x(2,npt)   = 0.0
  79.       x(2,npf+1) = 0.0
  80.       do i = 2,npf
  81.      j = npt+1-i
  82. c
  83. c        Calculate the camber line.
  84. c
  85.      xi = x(1,i)
  86.      xm = 1.0 - x(1,i)
  87.      if (xi.le.0.0 .or. xi.ge.1.0) then
  88.         y    = 0.0
  89.         dydx = 0.0
  90.      else
  91.         y    = (-ndata(1)/( 4.0*pi))*( xm*log(xm) + xi*log(xi) )
  92.         dydx = ( ndata(1)/( 4.0*pi))*(    log(xm) -    log(xi) )
  93.      endif
  94. c
  95. c        Calculate thickness.
  96. c
  97.      t = Nacat( xi,xb,yb,y2b,nbm,nb,ab,bb )
  98. c
  99. c        Finally, calculate node positions upon the foil surface.
  100. c
  101.      dy = sqrt(( t**2 )/( 1.0+( dydx**2 )))
  102.      dx = dy * dydx
  103. c
  104.      x(1,i) = x(1,i) + dx
  105.      x(1,j) = x(1,i) - 2.0*dx
  106.      x(2,i) = y - dy
  107.      x(2,j) = y + dy
  108.       enddo
  109. c
  110. c     Rescale dimensions if nessisary.
  111. c
  112.       call FoilScale( x,npt,scle )
  113. c
  114.       return
  115.       END
  116. c
  117. c-----------------------------------------------------------------------------
  118. c
  119.       SUBROUTINE Naca1b( naca,inaca,ndata )
  120. c
  121.       IMPLICIT none
  122.       INTEGER inaca
  123.       REAL ndata(3)
  124.       CHARACTER naca*(*)
  125. c
  126. c  This routine converts a NACA number into the parameters describing
  127. c  a foil.
  128. c
  129.       INTEGER icld,itau
  130.       INTEGER Ist
  131.       EXTERNAL Ist
  132. c
  133. c     Convert first four digits into integers.
  134. c
  135.       inaca= Ist( 2,naca(5:6),0 )
  136.       icld = Ist( 1,naca(4:4),0 )
  137.       itau = Ist( 2,naca(5:6),0 )
  138.       ndata(1) =   0.1*real(icld)
  139.       ndata(2) =  0.01*real(itau)
  140.       ndata(3) =  real(itau)/15.0
  141. c
  142.       return
  143.       END
  144. c
  145. c-----------------------------------------------------------------------------
  146. c
  147.       SUBROUTINE Naca6( x,npl,npmx,naca,inaca,scle,a,type,t_n )
  148. c
  149.       IMPLICIT none
  150.       INTEGER npl,npmx,inaca
  151.       REAL x(3,npmx),scle,a
  152.       CHARACTER naca*(*),type*3
  153.       LOGICAL t_n
  154. c
  155. c  This routine calculates the ordinates of a NACA 6 series foil.
  156. c
  157.       INTEGER nbmx
  158.       PARAMETER (nbmx=30)
  159. c
  160.       INTEGER npt,npf,nb,nbm
  161.       INTEGER i,j
  162.       REAL ndata(3),dx,dy,dydx,t,y
  163.       REAL ab,bb,rb,xb(nbmx),yb(nbmx),y2b(nbmx)
  164.       CHARACTER filnme*6
  165.       REAL Nacat
  166.       EXTERNAL Nacat
  167. c
  168.       npt  = npl+1
  169.       npf  = npl/2
  170. c
  171. c     Calculate the parameters describing the foil, from the given NACA
  172. c     number. Then calculate the cubic spline for the camber.
  173. c
  174.       call Naca6b( naca,inaca,ndata,filnme,type )
  175.       if (ndata(2).ne.0.0) then
  176.      call GetSpl(xb,yb,y2b,nb,nbm,nbmx,ab,bb,rb,filnme,t_n,ndata(3))
  177.       endif
  178.       if (.not. t_n) return
  179. c
  180. c     Set the x values according to a cosine distribution.
  181. c
  182.       call Cosdn( x,npl,npmx )
  183. c
  184. c     Calculate 'y' values for each value of 'x'.
  185. c
  186.       x(2,1)     = 0.0
  187.       x(2,npt)   = 0.0
  188.       x(2,npf+1) = 0.0
  189.       do i = 2,npf
  190.      j = npt+1-i
  191. c
  192. c        Calculate the camber line.
  193. c
  194.      call Naca6c( y,dydx,x(1,i),a,ndata(1) )
  195. c
  196. c        Calculate thickness.
  197. c
  198.      if (ndata(2).eq.0.0) then
  199.         t = 0.0
  200.      else
  201.         t = Nacat( x(1,i),xb,yb,y2b,nbm,nb,ab,bb )
  202.      endif
  203. c
  204. c        Finally, calculate node positions upon the foil surface.
  205. c
  206.      dy = sqrt(( t**2 )/( 1.0+( dydx**2 )))
  207.      dx = dy * dydx
  208. c
  209.      x(1,i) = x(1,i) + dx
  210.      x(1,j) = x(1,i) - 2.0*dx
  211.      x(2,i) = y - dy
  212.      x(2,j) = y + dy
  213.       enddo
  214. c
  215. c     Rescale dimensions if nessisary.
  216. c
  217.       call FoilScale( x,npt,scle )
  218. c
  219.       return
  220.       END
  221. c
  222. c-----------------------------------------------------------------------------
  223. c
  224.       SUBROUTINE Naca6b( naca,inaca,ndata,filnme,type )
  225. c
  226.       IMPLICIT none
  227.       INTEGER inaca
  228.       REAL ndata(3)
  229.       CHARACTER naca*(*),filnme*6,type*3
  230. c
  231. c  This routine converts a NACA number into the parameters describing
  232. c  a foil.
  233. c
  234.       INTEGER icld,itau,ibase,i
  235.       INTEGER Ist
  236.       EXTERNAL Ist
  237. c
  238. c
  239.       if (type(2:2).ne.'(') then
  240.      filnme = naca(1:3)//'0'//naca(5:6)
  241.      ibase  = Ist( 2,naca(5:6),0 )
  242.      i = -3
  243.       else
  244.      if (type(3:3).eq.'1') then
  245.         i = 0
  246.         filnme = naca(1:2)//type(1:1)//'00'//naca(4:4)
  247.      else if (type(3:3).eq.'2') then
  248.         i = 1
  249.         filnme = naca(1:2)//type(1:1)//'0'//naca(4:5)
  250.      endif
  251.      ibase  = Ist( 2,naca(4:4+i),0 )
  252.       endif
  253. c
  254.       inaca= Ist( 2,naca(8+i:9+i),0 )
  255.       icld = Ist( 1,naca(7+i:7+i),0 )
  256.       itau = Ist( 2,naca(8+i:9+i),0 )
  257.       ndata(1) =  0.1*real(icld)
  258.       ndata(2) = 0.01*real(itau)
  259.       ndata(3) = real(itau)/real(ibase)
  260. c
  261.       return
  262.       END
  263. c
  264. c-----------------------------------------------------------------------------
  265. c
  266.       SUBROUTINE Naca6c( y,dydx,x,a,cl )
  267. c
  268.       IMPLICIT none
  269.       REAL y,dydx,x,a,cl
  270. c
  271. c  Calculates the camber line for 6-series foils
  272. c
  273.       REAL pi
  274.       PARAMETER (pi=3.14159265358979)
  275.       REAL xm,am,g,h
  276.       REAL ax,da,ya
  277. c
  278.       xm = 1.0-x
  279.       am = 1.0-a
  280.       ax = a-x
  281.       if (x.le.0.0 .or. x.ge.1.0) then
  282.      y    = 0.0
  283.      dydx = 0.0
  284.       else
  285.      if (a.eq.1.0) then
  286.         y    = (-cl/( 4.0*pi))*( xm*log(xm) + x*log(x) )
  287.         dydx = ( cl/( 4.0*pi))*(    log(xm) -   log(x) )
  288.      else
  289.         if (a.eq.0.0) then
  290.            g = -0.25
  291.         else
  292.            g = (-1.0/am )*( a*a*( 0.5*log(a ) - 0.25 ) + 0.25 )
  293.         endif
  294.         h = ( 1.0/am )*( 0.5*am*am*log(am) - 0.25*am*am ) + g
  295.         if (x.ne.a) then
  296.            ya = ( 0.5*ax*ax*log(abs(ax)) - 0.5*xm*xm*log(xm)
  297.      &            +  0.25*xm*xm - 0.25*ax*ax )/am - x*log(x) + g - h*x
  298.            da = (    log(xm)*( xm + 0.5*xm*xm ) + 0.5*xm*xm
  299.      &            - log(abs(ax))*( ax + 0.5*ax*ax ) - 0.5*ax*ax
  300.      &            - 0.5*xm + 0.5*ax )/am - log(x) - h - 1.0
  301.         else
  302.            ya = ( -0.5*xm*xm*log(xm) + 0.25*xm*xm )/am
  303.      &            - x*log(x) + g - h*x
  304.            da = (    log(xm)*( xm + 0.5*xm*xm ) + 0.5*xm*xm
  305.      &            - 0.5*xm )/am - log(x) - h - 1.0
  306.         endif
  307.         y    = cl*ya/( 2.0*pi*( 1.0+a ) )
  308.         dydx = cl*da/( 2.0*pi*( 1.0+a ) )
  309.      endif
  310.       endif
  311. c
  312.       return
  313.       END
  314. c
  315. c-----------------------------------------------------------------------------
  316. c
  317.       SUBROUTINE GetSpl( xb,yb,y2b,nb,nbm,nbmx,ab,bb,rb,filnme,t_n,scl )
  318. c
  319.       IMPLICIT none
  320.       INTEGER nb,nbm,nbmx
  321.       REAL xb(nbmx),yb(nbmx),y2b(nbmx),ab,bb,rb,scl
  322.       CHARACTER filnme*(*)
  323.       LOGICAL t_n
  324. c
  325. c  Calculates the thickness spline for 1-series and 6-series foils
  326. c
  327.       INTEGER i
  328.       REAL y1,y2,xb2,xb3
  329. c
  330.       DATA  y2   / 3.0e35 /
  331. c
  332. c
  333.       call GetCor( xb,yb,nb,nbmx,rb,filnme,t_n )
  334.       if (t_n) then
  335.      nbm = nb-1
  336.      if (scl.ne.1.0) then
  337.         do i = 1,nb
  338.            yb(i) = yb(i)*scl
  339.         enddo
  340.      endif
  341.      xb2 = sqrt( xb(2) )
  342.      xb3 = sqrt( xb(3) )
  343.      ab  = ( yb(2)*xb(3)-yb(3)*xb(2) )/( xb(3)*xb2 - xb(2)*xb3 )
  344.      bb  = ( yb(2) - ab*xb2 )/xb(2)
  345.      y1  = bb + 0.5*ab/xb2
  346.      call Spln( xb(2),yb(2),nbm,y1,y2,y2b(2) )
  347.       endif
  348. c
  349.       return
  350.       END
  351. c
  352. c-----------------------------------------------------------------------------
  353. c
  354.       SUBROUTINE GetCor( xb,yb,nb,nbmx,rb,filnme,t_n )
  355. c
  356.       IMPLICIT none
  357.       INTEGER nbmx,nb
  358.       REAL xb(nbmx),yb(nbmx),rb
  359.       CHARACTER filnme*6
  360.       LOGICAL t_n
  361. c
  362. c  Reads in the coordinates of a NACA thickness distribution from file
  363. c
  364.       INCLUDE 'directory.h'
  365.       INTEGER ierr,i
  366.       REAL x,y
  367.       CHARACTER fullpath*128
  368.       LOGICAL fex
  369.       INTEGER Lnblnk
  370.       EXTERNAl Lnblnk
  371. c
  372.       fullpath = direct//filnme
  373.       inquire(file=fullpath,exist=fex)
  374.       if (.not.fex) then
  375.          print '(2a)', 'GetCor(): Cannot find file ',fullpath
  376.          t_n = .false.
  377.          return
  378.       endif
  379.       open( unit=10,file=fullpath,status='OLD',iostat=ierr )
  380.       if (ierr.ne.0) then
  381.          print '(2a)', 'GetCor(): Error opening file ',fullpath
  382.          t_n = .false.
  383.          return
  384.       endif
  385. c
  386.       read(10,*)
  387.       read(10,*)
  388.       x = 0.0
  389.       i = 0
  390.       do while( x.lt.100.0 )
  391.      read(10,*) x,y
  392.      i = i+1
  393.      xb(i) = 0.01*x
  394.      yb(i) = 0.01*y
  395.       enddo
  396.       read(10,*) rb
  397.       rb = 0.01*rb
  398.       nb = i
  399. c
  400.       close(unit=10)
  401. c
  402.       return
  403.       END
  404. c
  405. c-----------------------------------------------------------------------------
  406. c
  407.       FUNCTION Nacat( x,xb,yb,y2b,nbm,nb,ab,bb )
  408. c
  409.       IMPLICIT none
  410.       INTEGER nb,nbm
  411.       REAL x,xb(nb),yb(nb),y2b(nb),ab,bb
  412.       REAL Nacat
  413. c
  414. c  Returns an interpolated thickness for naca 1- and 6-series foils
  415. c
  416.       REAL y
  417. c
  418.       if (x.lt.xb(2)) then
  419.      y = ab*sqrt(x) + bb*x
  420.       else
  421.      call Spin( xb(2),yb(2),y2b(2),nbm,x,y )
  422.       endif
  423.       Nacat = y
  424. c
  425.       return
  426.       END
  427.